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^ — -  This  paper  contains  a  report  of continuing  investigat  100*: 
into  the  application  of  the  transition  state  theory  to  atom,  ion 
and  molecular  group  transfers. — Hi  Phis  papao  m  scacaaHata  our 
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sttMUda^on  the  formulation  of  Monte  Carlo  methods  of  simulation. 
We  find rh**-  _Tt  is  ^possible  to  derive  forms  for  the  rate  constant 
in  the  diabatlc  and  adiabatic  limits  which  lnmediateLy  suggest 
the  manner  in  which  Monte  Carlo  averages  should  be  carried  out. 


by  considering  the  reactive  inversion 
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We^i-Huat rata  the  methed 

of  a  bent  A-B-A  molecule  in  a  two-dimensional  system  of  disk- like 
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solvent.  -We  determine  the  activation  energy* for  the  transformation 


to  a  linear  transition  state,  and  v»  era "able  to  keep  a  tally  ed 
the  average  ediabaticity  of  the  reaction.  These  results  indicate 
that  the  extention  of  the  method  to  more  realistic  systems  ought 
to  yield  valuable  Information  about  the  mechanisms  of  reactions  in 
condensed  phases. 


X 
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lntroduct ion 

In  a  proceeding  paper*  the  transition  state  theory,  as 
formulated  by  Marcus,  was  examined  for  Che  purpose  of  finding 
algorithms  which  can  be  used  to  develop  computational  forms  for 
the  theory.  To  that  end,  an  arbitrary  system  was  viewed  as 
composed  of  a  definite  reactive  subsystem  which  is  surrounded 
by  an  environment.  The  local  microscopic  dynamics  of  the  reactive 
subsystem  can  be  examined  in  isolation.  Adjustments  in  these 
dynamics  can  be  accounted  for  as  the  local  subsystem  Is  placed 
into  the  remainder  of  the  solution.  Such  a  decomposition  of  an 
arbitrary  system  seems  to  imply  the  operation  of  weak  to  moderately 
weak  couplings.  However,  a  reactive  subsystem  can  be  defined  so 
that  a  sufficient  amount  of  solvent  is  included  with  the  reactant (s) 
to  encompass  all  species  which  experience  large  displacements  In 
the  course  of  a  reactive  transition.  Thus,  a  range  of  strengths 
of  interaction  between  the  reactant (s)  and  the  local  solvent  in  Che 
surroundings  can  bo  accounted  for. 

In  this  paper  we  continue  to  consider  the  construction  of  a 
computational  form  of  the  transition  state  theory.  We  examine  the 
transition  state  theory  In  a  form  which  makes  evident  the  mode  of 
application  of  the  Monte  Carlo  method  in  order  to  evaluate  certain 
average  values  such  as  the  energy  of  activation  and  the  pre~ 
exponential  factor.  By  examining  a  simple  two-dimensional  reactive 
system,  it  is  possible  to  see  the  emergence  of  some  of  the  limits 
on  and  limitations  of  a  Monte  Carlo  method  of  simulating  reactive 
systems 


3- 

Equl librium  Statistics  and  the  Development  of  the  Rate  Constant 

As  was  noted  in  Part  1,*  the  Marcus  transition  state  theory 
admits  Monte  Carlo  Techniques  for  the  purpose  of  carrying  out 
various  averages,  but  it  is  not  imnedlately  evident  how  this 
is  possible.  The  abstract  form  of  Marcus's  original  treatment 
disguises,  to  some  extent,  an  underlying  simplicity.  This 
simplicity  is  revealled  by  changing  to  a  more  model -dependent , 
less  abstract  formulation.  The  work  of  Glyde3  on  the  self- 
diffuslon  of  argon  suggests  the  essential  steps  one  needs  to  take 
in  order  to  establish  contact  with  the  Me  te  Carlo  methods  of 
averaging. 

In  the  following  paragraphs  the  transition  state  theory  is 
outlined  for  both  the  adiabatic  and  dlabatlc  limits.  This  sketch 
of  the  theory  lnnediately  suggests  the  proper  approach  to  take 
in  order  to  carry  out  the  Monte  Carlo  methods  of  averaging  or 
optimization. 

The  definition  of  an  adiabatic  transfer  was  considered  in 
Part  1.  If  a  reaction  proceeds  through  the  transition  state  via 
an  adiabatic  path,  the  transfer  species  then  occupies  a  state  of 
local,  transient,  and  atable  mechanical  equilibrium  for  the  Instate 
It  resides  In  the  configuration  of  the  transition  state.  Further¬ 
more,  an  adiabatic  path  for  a  reaction  is  one  on  which  the  transfer 
species  always  occupies  a  state  of  local,  transient,  and  stable 
mechanical  equilibrium  at  every  point  along  the  path.  In  other 
words.  In  the  adiabatic  limit,  fluctuations  in  the  conf iguratiosi 
of  the  entire  system  define  new  states  into  which  the  transfer 
species  moves.  This  movement  can  be  defined,  for  the  purpose 


f. 


-4- 


of  thermodynamic  argument,  to  take  place  reversibly.  Thus, 
the  t rans f ez -spec i es  adiabat ically  follows  a  migrating  position 
of  equilibrium.  Changes  in  the  configuration  of  the  environment 
carry  the  transfer-species  from  its  initial  to  its  final  location. 

The  diabatic  limit,  in  contrast,  applies  to  those  configura¬ 
tions  for  which  it  is  impossible  to  define  a  state  of  local, 
stable  mechanical  equilibrium  at  the  transition  point.  The 
migrating  species  must  then  tunnel  through  the  narrow  remaining 
barrier.  As  an  example,  one  can  consider  the  passage  of  an  atom 
or  ion  through  an  annular  opening  in  a  cryptate.4  If  the  effective 
radius  of  the  annulus  is  less  than  a  critical  value,  it  is  not 
possible  for  the  migrating  species  to  occupy  a  position  of 
equilibrium  at  the  centre  of  the  ring-like  structure.  However, 
if  the  effective  radius  increases,  a  stable  position  of  equilibrium 
can  result.  For  the  diabatic  limit,  therefore,  it  is  possible 
to  define  a  Bom-Oppenheimer  type  of  separation  of  the  local 
reactive  modes  from  the  remaining  modes  of  the  environment .  * 

As  will  become  evident  for  the  Monte  Carlo  simulations,  the  use 
of  a  Bom-Oppenheimer  type  of  separation  of  local  reactive  modes 
from  the  modes  of  the  environment  has  a  direct  bearing  only  on 
the  evaluation  of  Che  matrix  elements  which  are  associated  with 
the  transfer.  The  calculation  of  the  activation  energy  with  the 
use  of  Monte  Carlo  techniques  is  operationally  the  same  in  the 
diabatic  and  adiabatic  limits. 

The  total  Hamiltonian  operator  for  the  system  is  divided  into 
two  parts: 

«  *  Ht  +  Henv  (1) 


-S- 

for  which  H£  is  the  Hamiltonian  operator  for  the  transfer  species 
and  ^env  is  a8S°ciated  with  all  other  species.  The  operator  for 
the  transfer-species,  Ht ,  is  the  sum  of  the  kinetic  and  potential 
energy  operators . 

Ht  -  T  +  V.  ( 2> 

This  operator  is  one- dimensional  in  the  coordinate  which  coincides 
with  the  transfer-axis.  The  potential  energy  operator  V  spans 
the  Initial  and  final  configurations;  it  is,  for  example,  a 
potential  energy  function  with  a  double  minimum.  Solutions  in 
the  Bom-Oppenheimer  approximation  are  constructed  with  the  use 
of  basis  functions  and  *bfi  which  satisfy  the  local,  model 
Hamiltonian  operators 


(3) 


with 


Hb*b6  *  ‘b6*b4 


(4) 


and  the  indices  a  and  b  specify  the  physical  locations  of  the 
initial  and  final  states.  The  Indices  r  and  6  specify  energy 
states  for  the  transfer- species  at  the  locations  a  and  b. 

For  example,  and  can  be  local  harmonic  oscillator  potentials 
Thus,  and  are  the  usual  solutions  to  the  one- dimensional 


Schrbdinger  equation  for  the  oscillator.  In  terms  of  these  solutions 
and  in  terms  of  the  potential  energy  operators  V  and  V.  ,  one  has 


■  Hv  +  V  -  V. 


where  now  the  combinations  V  -  V  and  V  -  V.  operate  as  perturba- 


The  complete  eigenvalue  problem  is  written 


with  the  state  function  expanded  in  terms  of  a  basis  set  in  the 
»n<i  *h' 


In  the  usual  manner,  it  is  now  possible  to  find  sets  of  equations 
with  which  to  determine  the  expansion  coefficients  X^.  0°* 

finds 


H.  *  X,  -  [  L,  ,.X,. 

i>,iv  i y  iy.Ji  J* 


where  the  prime  on  the  smnation  excludes  i t  •  Jft.  The  L-matrix 
elements  are  defined  by 


L  -  y  c'1 

Lly.j6  k;tSi,,ktVkt>)5< 


here  S ^  is  an  element  of  the  inverse  overlap  matrix.  Roman 
indices  cover  the  initial  and  final  locations  (a  and  b)  while 
Creek  indices  cover  local  energy  states.  The  matrix  element  j f 


-  <kc|V  -  V.  |J6>  ♦  ‘k<IV  MS> 


In  defining  the  L-matrix  elements,  we  have  ignored  the  contribution 
from  the  kinetic  energy  operators  of  the  environment  which 
operate  in  the  space  of  the  transfer  species.  It  is  not  altogether 
clear  that  this  is  generally  permissible  for  vibrational  problems 
as  it  may  be  for  electron  transfer.5  It  is  of  course  relatively 
simple  to  amend  the  formalism  to  include  these  T-operators  in  the 
definition  of  the  L-matrix  elements,  indeed,  one  need  only  add 
7env  £o  V_VJ  et*n  l9*  •  the  diagonal  element  ^  is 

given  by 


Adiabatic  solutions  are  defined  by  the  solutions  to  the  equation 


H.  .  X.°  -  0. 

iy.iy  iy 


The  non-diagonal  terms  account  for  the  tunnel -transfer  of 

the  migrating  species  when  the  transfer  takes  place  In  the  diabetic 


From  the  first  order,  time- dependent  perturbation  theory,  the 
quantum  mechanical  transition  probability  Is  given  by 


.t 


W  "  W  il|f!eXP<‘tEi1)|Li>,f41;'i(Ei,  -  Ef«> 


(12) 


where  the  energies  E^  and  E^6  are  the  total  energies  of  the  entire 
system  in  its  initial  and  final  states,  and  y  and  6  specify  the 
state  of  the  transfer  species.  In  eqn  (12)  is  the  partition 
function 


l  exp( - sE . 
iy 


). 


(13) 


By  virtue  of  arguments  which  have  been  presented  in  detail  elsewhere,* 
it  is  possible  to  equate  the  transition  probability  in  the  diabatic 
limit  with  the  rate  constant: 


N  species  in  the  system.  is  now  a  classical  Hamilton  function 

for  the  system  in  its  initial  state.  The  energy  difference  in  the 
delta  function,  E^  -  *s  replaced  by  the  difference  in  the 

potential  energy  functions  in  the  initial  and  final  states  as  the 
kinetic  energy  terms  cancel.*'7 

The  integrations  over  the  velocities  in  Q  and  can  be 
carried  out  immediately.  As  there  is  no  Isolated  velocity 
in  Che  expression,  the  contributions  due  to  the  kinetic  energies 
cancel.  The  rate  constant  is 


■  r  h  «p<-'‘vi>m1.fJi2‘<vt  - 


Vf> 


(15) 


where  now  kd  stands  for  the  rate  constant  in  the  diabatic  limit.  The 
work  of  creating  the  initial  state  and  various  other  energies 
which  enter  the  activation  energy  appear  automatically  in  a  general 
analysis . 

We  now  assume  that  the  environment  (which  can  include  the  non- 
reaccive  vibrational  modes  of  the  molecular  framework  in  a  molecule 
which  undergoes  a  transformation)  can  be  treated  in  the  classical 
limit.  Therefore,  the  suxmnatlons  can  be  replaced  by  integrations 
with  the  result  that  the  rate  constant  now  is  expressed  as 


q  •  | dr  exp(-t>Vi) . 


(16) 


At  this  point,  we  make  the  assumption,  parallel  to  the 
assixnption  made  by  Clyde,*  that  the  major  contributions  to  the 
intagral  in  eqn  (16)  come  from  terms  in  the  exponent  which  are  near 
to  Vi0’  the  initial  equilibrium  potential  energy.  Thus,  on  expanding 
about  V, n.  one  finds 


*10 


+  *  K. 


(17) 


kd  *  fcjdr  dr  •*P<-»H1)|Ll>>f4l»«(V1  -  Vf> 


(14) 


(18) 


The  integration  involves  the  3N  velocities  and  coordinates  of  the 
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The  integration  over  r  yields  for  q 

q  -  (2./6)3N/2lAf*  exp(-6V10)  (19) 

where  jA|  is  the  determinant  of  the  matrix  of  second  derivatives. 
The  expression  for  the  rate  constant  is 

kd  -  ^<2,/e)3N/2|A|'  |  dr  expt-eCVj  -  V1Q)  1  |Lty  _  £4  |  2  4(Vt  -  Vf) 

(20) 

Write  the  delta  function  as^ 

‘<vi  -  V  '  TTfT4(xi.  '  x°f.)  <n> 

with 

l‘F|  -  |»(V1-Vf)/ax“| 0 j  (22) 


kd  ■  jrflTT<2’/t’>3N/2,A|4J‘is  »*P<-f(V(... -  V10)! 

‘23> 

where  ds  is  the  surface  element  in  the  reaction  hyperspace.  For 
this  analysis  it  Is  easily  defined;  it  la 

ds  -  lTdx“  <24) 

where  the  prime  excludes  the  single  coordinate  which  is  associated 
with  the  transfer  coordinate.  xj. 

We  assume  that  the  Condon  approximation  applies  The  matrix 
element  JL^  j4|2  now  can  be  removed  from  the  integration  Further, 
we  assume  that  the  major  contributions  to  the  Integrals  over  the 
x^  come  from  terma  in  the  vicinity  of  the  single  minimum  in  the 
transition  state.  Thus, 


which  is  the  difference  in  the  slopes  of  the  initial  and  final 
potential  energy  functions  evaluated  at  the  transition  point. 
Effectively,  this  quantity  is  the  difference  in  force  which  is 
experienced  by  the  transfer  species  es  it  passes  from  one 
potential  energy  surface  to  the  other  In  the  trenefer  region. 

The  single  coordinate  which  Is  Involved  in  the  Integration  over 
the  delta  function  is  x“.  The  Integration  which  Involves  the 
delta  function  has  the  effect  of  freezing  the  value  of  x“  for  the 
transfer  species  at  thet  value  which  corresponds  to  the  transfer 
point.  Thus. 


»(... «*..•>  * 


(2i) 


and  the  rate  constant  finally  is 


K  -  -T.-T- - exp{-e(V*  -  VlQ>). 


|4F| 


(26) 


By  way  of  contrast,  following  Clyde's  analysis  with  the 

modifications  Introduced  here,  one  finds  k  j.  the  rate  constant  in 

so 

Che  adiabatic  limit,  is  given  by 
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The  diabatic  mass  m^  is  just 
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(2»m  ) 


[A 


* 

3N 


expi - fi(VT 


rt)l- 


(27) 


It  is  necessary  to  note  the  fact  that  the  ratio  t  A  |  ^N/ |A*  |  j 
In  the  adiabatic  limit  will  not  be  the  same  as  the  similar  expression 
in  the  diabatic  limit,  cf.  eqn  (26).  However,  one  expects  the 
difference  to  be  (usually)  relatively  small. 

Beacuse  of  the  difference  in  dimensions  of  the  ratio  of 
generalized  force  constant  determinants  A.  it  is  easy  to  see  that 


|a|3N/|aMn-1  -  (28) 

is  effectively  a  single  force  constant.  The  frequency  is  an 
effective  frequency  which  can  be  assigned  to  the  transfer.  As 
Clyde  noted"*  for  the  self -diffusion  of  argon,  it  is  an  Einstein- 
like  frequency  With  the  restriction  that  the  ratio  of  force 
constant  determinants  is  essentially  the  same  for  the  diabatic 
and  adiabatic  limits,  it  is  possible  to  express  the  diabatic  limit 
as 


k,  -  (2  «m  ) 


T7T 


l*Tl 


—  «xp(-e<v+  -  v10)) 


3N-1 


-  (*  /md)  V£exp(-e(V^  .  V1Q) 


(29) 


The  effective  mass  is  then  defined  by 


2u/e\L 


it ,  f  * 


md  -  4me.  (31) 

An  Estimate  of  the  Pre-exponential  Factor 

The  simplest  estimate  of  the  size  of  the  pre-exponential  factor 
is  made  with  the  use  of  the  model  of  intersecting  harmonic  oscillator 

potential  energy  functions.  It  is  clear  from  our  other  investiga- 
1  8 

tions  '  that  the  diabatic  limit  occurs,  as  mentioned,  in  those 
instances  for  which  stable  mechanical  equilibrium  cannot  be  found  for 
the  transfer  species  in  the  transition  state.  It  is  equally 
clear  from  our  investigations**  that  these  conditions  depend  strongly 
upon  the  dominant  contribution  of  repulsive  forces  for  some  particular 
configuration  of  the  system.  Thus,  in  more  graphic,  physical 
terms,  it  is  easy  to  see  that  a  model  of  intersecting  harmonic 
oscillator  potential  energy  functions  is  an  adequate  representa¬ 
tion,  in  simple  terms,  of  an  actual  interplay  of  potential  energy 
functions  in  the  region  of  strong  repulsive  interaction 

For  the  simple  one -dimensional  oscillator,  the  potential  energy 
function  is 

V  -  %tm-?(z  -  to)'  (32) 

on  the  transfer  axis.  The  force  is 


(30) 


F  -  -  mw 2 ( z  -  z  0 ) 


(33) 


The  location  of  the  point  of  intersection,  a  cusp,  is  given  by 

*i»t  *  z*0  -  ^(a-Ocit.  <■-(-«•>  +  *'j%  <34> 


E D  is  the  difference  in  zero-point  energies  for  the  two  wells 
If  such  that  0*1,  then 


and  if  Et,  -  0,  then  zlnt  -  %Rab,  wh«re  Rab  is  the  distance  between 
the  locations  of  the  Initial  and  final  states. 

In  a  similar  manner,  one  finds  that  the  difference  F^  -  F ^  la 


Fj  -  Ff  -  +  zio  '  °ilf0) 


For  the  case  of  a  totally  symmetric  system, 


|  ftF  [  -  2«<frfL  c. 


where  c  is  and  to  is  the  distance  from  the  minimum  to 


We  found  in  Part  II  that 


In  terms  of  this  quantity,  and  in  tern  the  estimate  of  UF| 
for  the  intersecting  oscillator  potent  jne  finds  that  the 

effective  mass,  o>e.  Is  given  by 

2t2 

✓m  -  /m  -  —  .  (41) 

e 

In  view  of  the  definition  of  the  distance  t,  it  is  clear  from 
eqn  (41)  that  the  effective  mass  increases  dramatically  with 
increasing  distance  between  the  minimum  of  the  potential  and  the 
maximum  of  the  barrier.  Thi6.  of  course,  is  a  result  which  is 
well-established  in  the  lore  of  chemistry.  In  addition,  one  sees 
that  the  effective  mass  increases  with  temperature  as  tV  This 
result  can  be  attributed  to  the  increase  in  the  mean  displacement  of 
an  oscillating  particle  as  the  temperature  increases. 

Monte  Carlo  Methods  Implied  by  the  Formal  Expression  for  the  Rate 


To  begin,  a  reactive  transition  evolves  as  a  fluctuation  from 
the  thermodynamic  state  of  equilibrium  Thus,  one  prepares  the  syster 
in  an  Initial  atate  according  to  the  method  devised  originally 
by  Metropolis,  «t  el.9  Given  a  centrally  located  reactant  which 
is  surrounded  by  solvent,  the  entire  system  is  allowed  to  come  to 
equilibrium.  Each  atom  In  the  entire  system  of  reactant (s)  and 
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... .  lvcfit  individually  is  displaced  a  random  distance  along  an  axis. 

It  f.ltt  displacement  results  in  a  lower  overall  energy,  it  is  kept  and 
tt.c  sampling  proceeds  to  the  next  coordinate  and  to  subsequent  atoms. 
It,  on  the  other  hand,  the  displacement  results  in  an  increase  in 
the  energy,  the  displacement  is  kept  only  if  the  Boltzmann  weighting 
factor  is  less  than  a  randomly  generated  number  which  lies  in  the 

9 

range  of  zero  to  one  In  this  manner,  a  Markov  chain  is  constructed. 
The  process  of  sampling  continues  until  Che  average  energy  settles 
on  a  stable,  central  value.  At  this  point  the  system  is  in  an 
appropriate  initial  state  tor  the  reaction. 

Consider  as  an  example  a  simple  two-dimensional  reactive  system 
which  consists  of  a  non-linear  ABA  molecule  which  is  inmersed  in 
a  solvent  of  disk-like  atoms.  The  reaction  involves  a  simple 
inversion 

A  A 

h  •  B 
A  A 

in  much  the  same  manner  as  is  the  cast  for  the  inversion  of  ammonia.1^ 
The  preparation  of  the  initial  state  for  the  system  In  all  likely- 
hood  places  a  molecule  of  solvent  at  or  near  to  the  location  for 
tht  B-species  in  the  final  state.  Thus,  as  is  the  case  for  the 
solid  state.  here  as  well  it  is  necessary  to  form  a  vacancy  into 
which  b  can  migrate  The  work  which  is  required  to  form  this  vacancy 
is  part  of  the  energy  of  activation 

A  Monte  Carlo  simulation  of  the  formation  of  the  vacancy  can 
be  carried  out.  in  principle,  in  several  different  ways  The 


following,  however,  seems  to  be  the  simplest  and  HM»t  direct  Ur.c 
merely  continues  to  execute  Metropolis  samples  for  the  Monte  Carlo 
simulation.  The  simulation  continues  beginning  with  the  sv.tic  in 
its  initial  state  of  thermodynamic  (and  Monte  Carlo;  equilibrium. 

It  is  apparent  that  there  can  be  molecules  of  solvent  which  lie  in  the 
region  where  one  wants  to  create  a  vacancy.  A  test  therefore  is 
needed  to  see  if  any  molecule  must  be  moved  If  a  molecule  occupies 
the  space  where  a  vacancy  should  be,  then  it  is  necessary  to 
determine  whether  a  randomly  generated  displacement  will  move  that 
molecule  from  the  location  of  the  required  vacancy.  In  addition, 
the  move  can  only  occur  if  it  satisfies  the  criteria  imposed  for 
the  Metropolis  sampling  On  the  other  hand,  if  a  vacancy  already 
exists  or  nearly  exists  (£...,  a  molecule  lies  close  to  but  not  over 
the  location  of  the  vacancy),  then  it  is  necessary  to  test  each 
displacement  for  each  molecule  of  solvent  to  make  sure  that  no 
molecule  moves  into  the  location  of  the  vacancy. 

The  Monte  Carlo  simulation  is  carried  out  with  the  restriction 
that  a  vacancy  form  until  the  energy  of  the  system  stabilizes. 

The  state  for  which  there  exists  a  vacancy  intv'  which  the  transfer- 
species  can  migrate  therefore  defines  the  initial  state  for  the 
reactive  tranition. 

The  determination  of  the  energy  of  the  transition  state  is 
the  final  part  of  the  simulation  of  the  reaction.  Again,  there  art- 
several  ways  the  simulation  could  be  carried  out  One  wav. 
theoretically  the  soundest,  is  to  continue  to  generate  configurations 
for  which  the  reaction  is  allowed,  the  vacancy  remains 

For  such  a  conf igurat ion ,  assuming  a  linear  traiectorv  to  the 
final  state  through  the  transition  state,  the  t rans f er - spec l es  is 


placed  at  the  point  of  maximum  energy.  This  point  should  lie  between 
the  initial  and  final  positions  and  should  also  lie  within  the 
reactive  species  or  system  of  reactants.  With  the  reactive  system 
’'frozen"  in  this  transition-configuration,  the  environment  is  allowed 
to  relax  to  a  lowest  possible  energy.  The  energy  which  is  determined 
in  this  manner  corresponds  to  a  transition  state  energy  for  the 
single  configuration.  This  process  must  be  carried  out  a  sufficient 
number  of  times  in  order  to  get  a  statistical  sample.  Moreover, 
for  each  transition  state  configuration,  it  is  necessary  to  determine 
whether  the  transfer  species  is  in  a  state  of  mechanical  equilibrium. 
Such  a  calculation  clearly  consumes  much  more  computer  time  than 
either  of  the  two  preceeding  preparatory  computations. 

A  less  accurate,  approximate  approach  is  the  following,  its 
accuracy,  as  will  be  evident,  depends  upon  the  degree  to  which  the 
reactive  system  and  its  molecular  framework  is  insulated  or  isolated 
from  the  environment.  The  technique  is  simple.  Consider  the 
reactive  system  in  the  absence  of  the  surrounding  environment. 
Determine,  on  the  basis  of  reasonable  intuition  and  perhaps  an 
optimization,  the  location  of  a  probable  position  for  the  transfer- 
species  in  the  transition  state.  For  the  ABA  example,  a  reasonable 
choice  is  to  place  the  B-species  on  the  A-A  axis,  thus  creating  a 
linear  transition  configuration.  Now,  with  the  reactant  system 
"frozen"  in  this  configuration,  the  environment  is  allowed  to  inter¬ 
act  and  come  to  a  new  6tate  of  equilibrium.  The  energy  which  is 
found  in  this  manner  corresponds  to  the  energy  of  the  transition 
state . 

The  activation  energy  is  simply  the  energy  of  the  transition 
state  less  the  energy  of  the  original,  initial  state  of  equilibrium. 


In  spite  of  the  approximate  nature  of  this  approach  to  the 
determination  of  the  energy  of  the  transition  state,  it  is  possible 
that  some  adjustment  of  the  reactive  system  can  be  made  in  order 
to  obtain  a  more  accurate  value  of  the  energy.  For  the  ABA  system, 
for  example,  it  is  possible  that  a  slightly  bent  configuration 
corresponds  to  the  true  transition  state.  Such  a  calculation  can 
be  carried  out  and  tests  can  be  included  in  order  to  ensure  that 
the  Monte  Carlo  simulation  does  not  merely  regenerate  the  thermo¬ 
dynamic  initial  or  final  states. 

So  far,  only  a  brief  mention  has  been  made  of  the  determination 
of  mechanical  stability  for  the  transfer-species  in  the  transition 
state.  It  is  clear  that  for  the  last  approach  to  the  determination 
of  the  energy  of  the  transition  state,  each  conf igurat ion  which  is 
generated  as  the  system  stabilizes  about  an  energy  can  be  tested 
for  mechanical  stability.  Configurations  which  indicate  a 
mechanically  unstable  transition  state  correspond  to  cases  for  which 
a  diabatic  transition,  with  tunnelling,  applies.  It  is  a  simple 
matter  to  tally  the  diabatic  versus  adiabatic  configurations.  It 
is  a  less  simple  matter  to  remove  the  transfer-species  to  a  position 
on  the  reactive  trajectory  for  which  a  state  of  mechanical  equilibrium 
does  exist  end  from  there  to  determine  the  probability  for  tunnelling 
However,  such  calculations  can  be  carried  out.  It  is  possible, 
therefore,  to  obtain  an  appropriate  form  for  the  pre-exponential 
factor. 

One  should  note  the  fact  that  the  activation  energy  for  a 
configuration  for  which  a  diabatic  transfer  holds  must  be  determined 
still  with  the  transfer  species  In  the  transition  state  conf igurat ion. 
even  though  this  configuration  is  mechanically  unstable. 
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An  interesting  result  of  these  simulations  is  the  statistical 
definition  of  an  adiabatic  ;>«.  rv  ue  diabatic  limit  for  a  reaction. 

Consider  an  individual  configuration  for  the  transition  state 
From  that  which  has  been  stated  here,  and  in  part  I,  it  is  clear  that 
the  transfer  corresponds  to  one  either  in  the  adiabatic  or  diabatic 
limit.  At  the  microscopic  level,  there  is  no  gradation  of  diabatlcity, 
the  reaction  is  either  diabatic  or  It  is  adiabatic.  However, 
in  a  statistical  sense,  for  a  given  reactive  system  it  is  possible 
that  the  Monte  Carlo  simulations  will  generate  configurations  for 
which  a  fully  adiabatic  transfer  applies  part  of  the  time  and  for 
the  remainder  of  the  conf lguratlons  a  diabatic  limit  applies.  The 
overall  value  of  the  pre-exponential  factor  therefore  must  reflect 
a  statistical  balance  between  the  two  extreme  limits.  Thus,  one 
sees  that  a  given  reactive  system  can  exhibit  a  type  of  transfer 
which  corresponds  overall  to  a  statistical  weight  of  the  diabatic 
against  the  adiabatic  limit.  That  this  is  indeed  the  case  is 
already  evident  in  a  Landau-Zener  type  of  analysis. The 
difficulty  in  applying  the  Landau-Zener  analysis  to  any  system 
which  undergoes  a  reaction  which  belongs  between  the  extremes  of 
diabaticity  and  adiabaticity  is  well  known.  Monte  Carlo  methods  of 
simulation  and  averaging,  on  the  other  hand,  provide  the  means  to 
investigate  the  shades  of  diabatlcity  which  one  expects  for  real 
systems 

Algorithms  for  Determining  Mechanical  Stability  in  the  Transition 

State  ) 

One  can  test  for  mechanical  stability  of  the  transfer-species 


in  the  transition  state  configuration,  the  test  is  carried  out  in 

the  classical  limit.  It  is  necessary  simply  to  determine  the 

sign  of  the  second  order  coefficient  in  the  Taylor  expansion 

for  the  motion  of  the  particle  along  the  transfer  axis.  The  expansion 

uses  pair-wise  interactions  between  the  transfer-species  and  each 

atom  of  the  reactant (s)  and  the  solvent. 

The  symaetry-adaptable  Taylor  aeries  for  a  general  scalar 
12 

function  g(r)  Is 

g(r+R)  *  •'It  J  (rn/nl)[(-i)ntlA  P  (R.r)I  (R)  (42) 

'  -  n-0  t  nt 

in  which*^ 

An|  «  0  for  l  >  n  and  n  -  t  odd 

*  for  «  i  n  and  n  -  >  even  <4J) 

and 

In,(R)  -  —i—l  dk  kn+2f(k)j  (kR)  (44) 

nl  (2i) *' o 

In  eqn  (44)  j^fkJl)  is  the  spherical  Bessel  function  of  the  first 
kind.13  and  f(k>  ia  given  by12 

f(k)  ■  4* |  dr  r2g(r)J D(kr) .  (45) 

*  o 

In  eqn  (42).  Pt(R*r)  is  the  Legendre  polynomial  of  order  i  and 

R*r  is  the  scalar  product  of  the  unit  vectors  for  R  and  r.  it  is  the 
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the  cosine  of  the  angle  between  the  vectors. 

When  the  z-axis  is  specified  to  be  the  transfer  axis,  eqn  (42) 
reduces  to 


g(r+R)  •  l  l  Jr(-i)n 

n-0  nI 


AntPi<cosVInt<R> 


(46) 


The  general  expansion  of  a  potential  of  the  form  1/r"  is 


—To  lN;1(«-.1)<-/')n 


(50) 

The  condition  for  stability  is  derived  simply  from  eqn  (50).  it  is 


where  now,  dr  is  the  angle  between  the  location  of  the  source  and 
the  z-axis. 

Mechanical  stability  for  the  motion  along  the  z-axis  at  the 
location  of  the  transition  state  is  guaranteed  if 


[5  +  16P. (coseR^)  j  At  -  (22  +  56P,  (cos*^) 


(51) 


g<(R)  -  (-i)t+2A2tPt(coseR)I28(R) 

-  -  ^  (1J0<R>  ■  »2<CO»eR)l?2<R>j 


<47) 


An  expression  of  this  type  is  no  more  difficult  to  evaluate,  and 
consumes  about  the  same  amount  of  computer  time,  as  the  evaluation 
of  each  of  the  individual  pair-wise  contributions  to  the  energy 
of  the  system. 


is  greater  than  zero. 

In  order  to  determine  if  the  transfer-species  occupies  a  state 
which  is  mechanically  stable  ir.  the  transition  state,  it  is  necessary 
to  evaluate  g2(RA)  for  each  species  in  the  system.  Thus, 

g,  *  <(R,)  >  0  (48) 

i  ’  J 

where,  clearly,  the  prime  on  the  sunmacion  excludes  the  transfer- 
species.  As  an  example,  consider  a  system  of  particles  associated 
through  the  interactir,.  of  the  Mie  potential: 

J 

(49) 


Computational  Aspects  of  the  Monte  Carlo  Simulation  of  a  Reactive 
Transfer  in  the  ABA  System 

The  system  is  two-dimensional  and  consists  of  a  bent  A-B-A 
molecule  which  la  immersed  in  a  solvent  of  uniform  disks  The 
origin  of  the  coordinate  system  Is  located  initially  at  the  mid¬ 
point  of  the  A-A  axis.  The  B-specles  is  located  at  a  perpendicularly 
displaced  point  over  this  midpoint.  The  solvent  is  distributed 
about  the  solute  in  concentric  circles  the  radii  of  which  are 
multiples  of  the  solvent  diameter;  the  first  radius  includes 
an  effective  radius  for  the  solute.  No  form  of  optimization  is  used 
to  construct  this  initial  configuration. 


V(r)  -  -  A/  r*  +  B/r>‘ 
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The  system  is  equilibrated  first  at  4K ,  this  step  is  suggested 

14 

in  a  similar  treatment  which  was  carried  out  by  Simons  The 

equilibration  is  accomplished  in  approximately  60.000  cycles  of 
the  Monte  Carlo  simulation.  The  temperature  is  then  raised  to 
300K  The  simulation  is  allowed  to  proceed  through  approximately 
20.000  cycles.  At  that  point,  memory  of  the  initial  energy  and 
configuration  of  the  system  is  erased  and  a  new  equilibration  la 
sought  in  approximately  60,000  cycles.  At  the  end  of  the  equilibra¬ 
tion  at  300K,  the  system  is  considered  to  be  prepared  in  an  Initial 
state.  It  Is  evident  that  the  system  has  reached  a  state  of 
equilibrium  at  a  given  temperature  when  the  fluctuations  in  the 
energy  diminish  and  an  essentially  stable,  almost  constant  value 
of  the  energy  is  reproduced  every  thousand  cycles  or  so. 

Although  emphasis  has  been  placed  on  the  process  of  forming 
a  vacancy  into  which  the  migrating  species  can  transfer,  this  Is 
not  altogether  a  necessary  step  in  the  simulation.  It  is  possible, 
in  fact,  merely  to  seek  good  configurations  for  the  transition 
state.  The  activation  energy  is  simply  the  difference  between  the 
absolute  energy  of  the  transition  state  and  the  energy  of  the 
initial  state.  Nevertheless,  it  appears  to  hasten  the  convergence 
of  the  simulation  to  optimal  configurations  and  values  for  the 
energy  if  the  step  to  find  the  vacancy  is  included.  This  process 
is  carried  out  as  follows.  To  begin,  a  test  is  carried  out  to  see 
if  a  molecule  of  solvent  occupies  a  position  in  the  eventual  location 
of  the  migrating  B-species.  If  the  test  proves  to  be  positive,  the 
molecule  of  solvent  is  displaced  along  the  y-axis  (the  transfer  1 

direction)  until  an  appropriate  void  is  created.  With  the  molecule 
of  solvent  held  rigidly  in  this  position,  the  remaining  solvent  is 
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allowed  to  equilibrate  about  the  solute  and  its  vacancv  wiih  the 
temperature  at  300K 

Once  the  vacancy  is  created,  an  Initial  attempt  to  form  the 
transition  state  is  made.  The  transition  state  is  prepared  simply 
by  moving  the  B-species  to  the  origin  of  coordinates  at  the  midpoint 
on  the  A-A  axis  The  A-species  are  allowed  to  equilibrate  only 
along  the  bond  axis,  this  is  done  in  order  to  prevent  the  re-establish¬ 
ment  of  the  Initial  state  or  the  establishment  of  the  final  state 
The  solvent  equilibrates  as  usual.  The  entire  process  is  again 
carried  out  at  300K.  With  the  generation  of  each  configuration  in  the 
process  of  equilibration,  the  transition  state  is  tested  to  see  if 
the  transfer  species  occupies  a  position  which  Is  mechanically 
stable  If  the  transfer  species  alone  occupies  a  position  within 
a  configuration  and  the  species  is  mechanically  stable  in  that 
position,  the  transfer  is  termed  adiabatic.  The  average  adiabaticity 
is  determined. 

The  running  average  of  any  quantity  is  kept  with  the  use  of 
the  simple  recursion  relation: 

rN  -  »(fN  +  (N-l)V,)  (52) 

where  f^  is  the  value  of  F  (any  measurable  quantity)  for  the  N-th 
individual  configuration,  and  Fjj  j  and  F^  are  the  average  values 
after  N-l  and  N  cycles  respectively.  With  the  use  of  this  recurrence 
relation  for  the  average,  it  is  easy  to  see  the  emergence  of  a 
relatively  stable  value  for  any  quantity  which  is  averaged  It 
is  also  easy  to  see  the  effect  of  long-term  memory,  especially  when 
the  simulation  is  programed  to  evaluate  averages  over  smaller 
aubcol lections  of  randomly  generated  configurations 


A 
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In  addition  to  the  details  of  the  calculation  which  have  been 
described  above,  we  included  a  circular  boundary  with  a  radius  equal 
to  the  distance  of  the  outermost  molecule  of  solvent  from  the 
coordinate  origin  plus  one  solvent  diameter.  The  interaction 
between  any  species  and  the  wall  of  this  container  was  taken  to 
be  infinitely  repulsive. 

The  potential  energy  which  was  used  throughout  the  simulation 
was  the  Lennard-Jones/Mie  potential: 

V(r)  -  D!<r0/r)'-  -  2(rt,/r)‘l  (53) 

where  D  is  the  dissociation  energy  and  r0  is  a  diatomic,  equilibrium 

bond  distance.  The  dissociation  energy  for  the  A-B  bond  was 
- 19 

taken  to  be  5.0-10  J  and  the  dissociation  energy  for  the 

-19 

A-A  was  half  that  value,  2.5*10  J.  For  simplicity,  the  dissociation 

energy  for  the  interaction  between  solvent  and  between  each  atom 
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of  the  solute.  A  and  B.  and  the  solvent  was  taken  to  be  9*10  J. 

The  equilibrium  distance.  r0,  for  the  A-B  bond  was  kept  a  constant 
O.lnm.  In  contrast,  the  bond  distance  for  the  two  A-specles 
was  allowed  to  run  from  a  minimum  value  of  O.lnm  (allowing,  in 
vacuum,  an  almost  equilaterlally  triangular  molecule)  to  0.2nm. 

The  solvent- solvent  and  solvent-solute  atom  distance  was  taken  to 
be  0.3nra  which  is  approximately  the  van  der  Waals  separation  for 
water. 

Results  of  the  Calculation 


energy  decreases  as  the  distance  between  the  two  A-specles  increases, 
it  becomes  energetically  less  costly  to  move  the  migrating  h-species 
into  the  configuration  of  the  transition  state  As  the  data  indicate 
when  the  A-A  bond  distance  is  allowed  to  equilibrate  in  the  transi¬ 
tion  state,  the  result  is  generally  an  adiabatic  transfer  It  is 
clear  that  the  adlabaticity  of  the  transfer  increases  as  the 
fundamental  A-A  distance  increases. 

The  question  naturally  arises  with  these  types  of  simulations 
as  to  whether  the  number  of  molecules  of  solvent  which  surround 
the  solute  is  sufficiently  large  to  account  accurately  for  the 
effect  of  the  solvent  on  the  activation  In  order  to  test  the 
adequacy  of  the  account  of  the  solvent,  an  additional  layer  was 
considered.  The  result  of  the  simulation,  which  took  more  than 
twice  the  time  on  the  computer  the  original  calculation  cook,  was 
identical  to  the  result  with  the  smaller  sample  of  solvent.  It 
appears,  for  the  size  of  solute,  A-B-A.  which  we  considered,  two 
distinct  layers  of  solvent  in  concentric  circles  are  sufficient 
to  estimate  the  energy  of  activation. 

The  calculations  also  reveal  more  than  Table  1  alone  shows 
It  1 s  clear  in  examining  the  evolution  of  the  adiabaticity  as  the 
system  moves  toward  an  optimal  value  for  the  transition  state 
that  the  initial  adiabaticity  is  small;  this  is  illustrated  in 
Table  2.  Adiabaticity  grows  as  the  A-A  bondiength  expands  to 
accomodate  the  enclosed  B-specles.  This  occurence  suggests  that 
if,  for  the  reason  of  a  constraint  due  to  a  larger  and  more 
)  complicated,  effectively  rigid  molecular  geometry,  a  molecule  cannot 

easily  expand  to  accomodate  a  transfer  species  in  a  state  of 
mechanical  equilibrium  in  the  configuration  of  the  transition 


As  intuition  and  experience  leads  one  to  suspect,  the  activation 
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j, t a 1 1* .  then  the  heavy-centre  transfer  will  tend  to  be  diabatic. 

The  relative  inflexibility  of  a  molecular  framework,  therefore, 
is  surely  an  obvious  limitation  on  the  adiabaticity  of  any  heavy- 
centre  atom  or  molecular  group  transfer. 

Discussion 

It  is  appropriate  to  consider  at  least  briefly  the  comparison 
between  the  theory  which  is  developed  in  this  paper  and  the 
transition  state  theory  of  reactions  in  the  gas  phase  as  formulated 
by  several  people,  especially.  Miller*^  and  Pechukas.*6 

First  of  all,  an  assumption  is  made  here  that  there  is  sep¬ 
arability  in  the  quant al  limit  in  the  sense  of  the  Bom-Oppenheioer 
approximation.  This  is  an  essential  assxanption,  for  it  prescribes 
the  form  which  the  expression  for  the  transition  probability 
(and  hence  the  rate  constant)  must  take.  The  development  of 
each  expression  for  the  rate  constant  in  the  diabatic  and  adiabatic 
limits  proceeds  strictly  governed  by  the  limitations  which  are 
dictated  by  assuming  the  classical  limit  for  all  degrees  of  freedom 
which  are  associated  with  the  solvent.  The  single  quantum  mechanical 
contribution  to  the  expression  for  the  rate  constant  in  the  diabatic 
limit  arises  from  the  need  to  consider  tunnelling  across  t he 
residual  barrier  to  the  transfer. 

In  principle,  all  the  difficulties  which  are  Inherent  in  the 
gas  phase  theory*6  of  the  transition  state  are  also  contained  in 
the  theory  for  the  condensed  phases  as  well.  Because  we  have 
chosen  to  examine  the  classical  limit  for  the  environmental 
contributions ,  the  expression  for  the  rate  constant  in  that  limit 


is  well-defined,  this  is  also  the  case  for  the  gas  phase.*6  As 
long  as  the  Bom-Oppenhelmer  approximation  applies  to  a  transfer, 
the  transition  probability  (from  second  order  time-dependent 
perturbation  theory)  also  applies.  The  classical  limit  for  the 
environmental  contributions  also  applies.  However,  in  the  limit 
of  strong  interactions,  where  quantal  limits  al6o  apply,  the 
theory  which  is  developed  here  is  inappropriate  Th  ,s ,  it  is 
Incorrect  to  attempt  to  apply  our  results  to  highlv  endo-  or 
exothermic  reactions  It  is  well-known,  also,  from  the  theory  of 
the  electron  transfer  reaction  that  a  more  appropriate  (but .  perhaps 
not  yet  entirely  accurate)  form  for  the  rate  constant  arises  from 
the  proper  consideration  of  the  saddle  point  Integration  which 
enters  the  theory.*^  Such  an  analysis  has  not  yet  been  carried 
out  for  these  heavy-centre  transfers  Moreover,  it  is  not  iasnediatel y 
clear  how  such  an  analysis  should  be  carried  out  in  order  to  derive 
expressions  for  the  rate  constant  which  lend  themselves  inmediately 
to  Monte  Carlo  simulation. 

It  is  evident,  therefore,  that  the  results  which  we  have  presented 
for  the  simulation  of  rate  constants  for  heavy-centre  transfers 
In  the  condensed  phase  must  be  considered  in  the  context  of  the 
limits  which  define  the  problem  The  transfers  which  we  are  able 
to  consider  are  those  for  which  the  enthalpy  of  reaction  is  in  the 
so-called  "normal”  range;  the  reactions  can  be  neither  highly 
exothermic  nor  highly  endothermic.  Finally,  it  is  mandatory  that 
the  reactions  be  considered  to  take  place  in  the  classical  limit 
with  reference  to  all  the  possible  dynamical  contributions  from  the 
solvent . 
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Table  1 

A.  (2.5/0.100/5.0/0.100)* 

E0  -  -54.87  (at  300K) 

E  -  -53.85 
v 

Et  -  -5165 

E  -  3.22  (448  kJ  nol-1) 

a 

(final  adiabaticity .  0.98) 

B.  (2.5/0.125/5.0/0.100) 

E„  -  -54  36 

Ev  -  -53  82 
Et  -  -51.55 

E4  -  2  86  (439  kJ  mol'1) 

( final  adiabaticity:  0.88) 

C  (2. 5/0. 180/S. 0/0. 100) 

Ed  -  -54.42 
Ev  -  -53.85 
Et  -  -  53.43 

Ea  -  0.99  (153  kJ  mol'1) 
(final  adiabaticity:  0.88) 
D  (2.5/0.200/5.0/0.100) 

Eo  -  -54.33 
Ev  -  -54.27 
Ec  -  -54.20 

Ea  -  0.13  (20  kJ  nol’1) 
(final  adiabaticity:  0.99) 


♦The  notation  (daA/raA/DAB/RAB*  lndlcat«9  the  p«r»®eters  used  for 


the  solute:  is  the  A-A  dissociation  energy.  is  the  A-B 

bond  dissociation  energy.  is  the  A-A  equilibritsc  pairwise 

separation  In  the  gas  phase,  and  finally  R^g  is  the  equi librium 
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A-B  separation.  The  energies  are  given  in  Joules*10  and  the 
distances  in  nm. 

Table  2  The  growth  of  adiabaticity  for  (2.5/0.180/5.0/0/100) 
Cycles  E  Adiabaticity 


2,625 

-52.27 

0.0015 

4.375 

-52.89 

0.  39 

6.125 

-53.12 

0.48 

7.875 

-53.22 

0.59 

9.625 

-53.33 

0  66 

11.375 

-53.36 

0.71 

13.125 

-53.38 

0.75 

14.875 

-53.38 

0.78 

16.625 

-53.40 

0.80 

18.375 

-53.43 

0.82 
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